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Abstract. Broadband emission from relativistic outflows (jets) of active galactic nuclei (AGN) 
and gamma-ray bursts (GRBs) contains valuable information about the nature of the jet itself, 
and about the central engine which launches it. Using special relativistic hydrodynamics and 
magnetohydronamics simulations we study the dynamics of the jet and its interaction with the 
surrounding medium. The observational signature of the simulated jets is computed using a 
radiative transfer code developed specifically for the purpose of computing multi-wavelength, 
time-dependent, non-thermal emission from astrophysical plasmas. We present results of a series 
of long-term projects devoted to understanding the dynamics and emission of jets in parsec-scale 
AGN jets, blazars and the afterglow phase of the GRBs. 



1. Introduction 

Relativistic jets are found in a variety of astrophysical scenarios such as active galactic nuclei 
(AGN), microquasars, gamma-ray bursts (GRBs) and those resulting from tidal-disruption 
events (TDEs). While jet temporal and spatial timescales can vary from relatively small and 
short-lived (GRB jets) to very large and long-lived (AGN jets), what they have in common is that 
they are relativistic, seem to be very collimated and appear to be launched from a region very 
close to the accreting black hole (see e.g., [1] for a recent review and introduction to relativistic 
jets). 

In this paper we review a series of projects whose aim is to understand the dynamics and 
emission from relativistic jets by focusing on specific scenarios appearing in relativistic jets 



of AGNs, blazars and GRBs. We begin by describing the numerical method in section [2j 
where we describe the equations of special relativistic magnetohydrodynamics and the code 
MRGENESIS (section 12. ip . and then give some details about the treatment of non-thermal 
particles (section l2.2p and the radiative transfer (section l2.3p . both of which are handled by the 
SPEV. In section [3] we present three of the applications of MRGENESIS+SPEV method, and 
we conclude in the section [H 



2. Numerical method 

The numerical method we use consists of three parts, each dealing with a different aspect of the 
jet dynamics and emission. In the section I2TT1 we discuss the equations of the special relativistic 
magnetohydrodynamics (RMHD). In sections 12.21 and 12.31 the treatment of non-thermal particles 
and emission, as well as that of the radiative transport is explained. Finally, the structure of 
the SPEV code is explained in 12.3.11 



2.1. Special relativistic magnetohydrodynamics 

We solve the system of equations of ideal special RMHD [2], which consists of the mass and the 
total-energy momentum conservation^] 

V a (pu a ) = ; V a T a P = Q, (1) 

where the covariant derivative is denoted by V Q . p is the fluid rest-mass, u a is the fluid 4- velocity, 
and the energy-momentum tensor is defined as 

T aP := ph*u a U P + V *ri aP - b a b? . (2) 

The Minkowski metric is rf 1 ^ = diag{—\, +1, +1, +1), and the total pressure isp*:=p+|6| 2 /2, 
where p is the fluid pressure. b a is the magnetic field four- vector, 

b° := T (v • B) , b i :=^ + v% Q . (3) 

Here B is the magnetic field 3-vector and v the fluid velocity 3-vector, whose Lorentz factor is 
defined as T := (1 — v 2 ) 1 ^ 2 . The total specific enthalpy is defined as 

h*:=h+^ = l + e + V + ^, (4) 
P P P 

where e is the specific internal energy, which depends on the equation of state p = p(e, p). We use 
the TM approximation to the Synge equation of state, where the specific enthalpy h is defined 
as [3] h = (5/2)(p/p) + \J (9/4)(p/p) 2 + 1. Finally, the magnetic field has to satisfy the induction 
equation 

9B 

— -Vx(vxB) = 0, (5) 

and the divergence constraint V • B = 0. 

We use the code MRGENESIS [U El El [7] to numerically solve the system of equations 
of RMHD. MRGENESIS uses operator-splitting, finite- volumes, explicit method and employs 
approximate Riemann solvers for the computation of intercell fluxes, as well as total variation 
diminishing second and third-order Runge-Kutta methods for the time integration (for a recent 
overview of jet simulation codes see e.g. [8j). MRGENESIS has been massively parallelized 
using OpenMP and MPI libraries and has achieved reasonable scaling up to 10 4 cores on the 
MareNostrum machine at the Barcelona Supercomputing Centre. 



1 In this section we assume a system of units where the speed of light c — 1 and the factor l/\/47r is included in 
the definition of the magnetic field. 



2.2. Non-thermal particles and emission 

As discussed in the section [TJ a small fraction of very energetic particles ("non-thermal electrons" 
in the following) is responsible for the observed jet emission. In [9] we discuss in detail the 
method which deals with the spatial and temporal evolution of the non-thermal electrons in a 
magnetized relativistic fluid. Here we summarize the most important properties of the method. 



2.2.1. Non-thermal electron evolution We assume that electrons are injected at relativistic 
shocks and use the phenomenological injection term which denotes the number of particles 
injected per unit volume and unit time, 

Q(7) = Qo (—^—] S (7; 7min,7max) , (6) 
V7min/ 

where 7 is the electron Lorentz factor 7 = E/(m e c 2 ), E being its energy and m e being its mass. 
Qo is the normalization of the injection energy spectrum, q is the power-law index and the 
function S(x; a, b) has the value 1 if a < x < b and otherwise!. The temporal evolution of the 
electron energy distribution in the fluid comoving fram^f] is governed by the kinetic equation 
[TO]: 

^|^ + ^[7n(t;7)] = 0(7), (7) 

where n(t; 7) is the number density of electrons in an interval d7 around the Lorentz factor 7 
at the time t. Energy gains and losses are described by the term 7 := d7/dt. For sufficiently 
small intervals of time we can write the energy losses as 

^L = k al -k s7 2 , (8) 

where the adiabatic gains or losses are described by the quantity 

dlnp 



dt 



which denotes the compression or expansion of the fluid element, while the synchrotron losses 
are denoted by 

Aa T B' 2 



3m 2 . c 2 



where B' is the magnetic field strength in the fluid comoving frame and o~t is the Thomson cross 
section. 

As described in section 3.2 of [S], we discretize the electron energy space in a number of 
logarithmically spaced bins and use a semi-analytic solver to accurately solve the equation ([7]) 
with a minimum computational effort. 



2.2.2. Non-thermal radiation processes We consider the following radiation processes in our 
simulations: synchrotron emission, synchrotron-self Compton (SSC) scattering, and the external 
inverse-Compton (EIC) scattering. 

Synchrotron emission is produced when high-energy electrons gyrate around the magnetic 
field lines, producing a non-thermal, broadband emission spectrum. It is computed for each 

2 We fix q = 2.3 and determine the rest of the parameters by assuming a proportionality between the number 
and energy density of the non-thermal electrons and those of the thermal fluid (see equations 36 and 37 of [5]). 

3 We assume that non-thermal electrons are advected with the thermal fluid. The spatial distribution of non- 
thermal electrons is represented by Lagrangian tracer particles. 



energy bin of the electron distribution using the interpolation method described in the sections 
2.1.3 and 4.3.1 of [II]. 

Inverse-Compton scattering is a scattering of a low-frequency photon off a high-energy 
electron, whereby the frequency of the outgoing photon can be many orders of magnitude larger 
than that of the incoming one. In our numerical method it is computed analytically (see section 
2.2.2 of [H]) assuming that the incoming emission spectrum is a power-law, as is the energy 
spectrum of the electrons off which the photons scatter. This method has been used in [T2] to 
compute the EIC scattering of the ultraviolet stellar, photons, as well as in [T3] to compute both 
the SSC and EIC emission in the blazar jets. 

2.3. Radiative transfer 

To compute the time- and frequency-dependent synthetic image we use the algorithm described 
in the Appendix A of [9] . Both this algorithm and the one described in the previous section are 
building blocks of the code SPEV. SPEV solves the radiative transfer equation 



where I u , j u and a u are the specific intensity, emissivity and absorption coefficient at frequency 
v. s is the distance along the line of sight which is seen by the observer at the time T. The 
relation between the time of observation and the time t in the laboratory frame (the frame in 
which the jet evolution is simulated) is given by 



where D is the distance of the jet injection point from the Earth and we choose that s = at 
that point. As can be seen, in order to solve the equation (j^J), it is not enough to process each 
saved simulation snapshot at time Snapshot- Instead, it is necessary to process all snapshots 
simultaneously because of the fact that distance from the observer and the time at which the 
snapshot has been taken are not independent, but must satisfy the equation (|10p . In practice, 
this means that the radiative transfer equation ([9|) cannot be solved before all the information 
along the line of sight has been gathered. This makes the problem tightly coupled and highly 
non-local, complicating the code paralellization and increasing the CPU and memory costs. 

2.3.1. Structure and paralellization of the SPEV code The evolution and emission from the non- 
thermal particles (section 12. 2[) and the computation of the resulting emission (section I2.3P are 
done by the SPEV code whose structure can be sketched as follows: after creating the initial 
Lagrangian particles (representing the pre-existing non-thermal electrons), SPEV iteratively 
runs through the three phases: evolution of existing particles, injection of new particles, and 
calculation of the emission (i.e., contribution to the pixels of a virtual detector on which the 
emission is computed). 

SPEV has two important characteristics very useful for its paralellization: the evolution of 
each Lagrangian particle during the simulation does not affect the evolution of the rest of the 
particles, and the contribution to each pixel of the detector is independent from the contribution 
to any other pixel. Therefore, SPEV is parallelized over particles and detector pixels. The 
data is partitioned into different sets and distributed among different computing units (e.g., 
multiple multicore-nodes connected by means of the network or by means of buses). Each set 
of particles can evolve separately without the need of any kind of communications during the 
iterative process. When the iterative process ends, data is transferred to the master node and 
then the radiative transfer equation ([9]) is solved for each pixel. More details about the parallel 
implementation are provided in [14^ 115]. This parallel approach has shown very good scalability 
on different HPC-systems, especially on medium scale multi-sockets multicores of up to 50 cores. 




0) 



t = T + (D- s)/c 



(10) 



3. Applications 

In this section we outline a number of application of MRGENESIS + SPEV approach to 
relativistic jets. In section 13.11 we discuss the radio emission from parsec-scale jets. The 
section 13.21 is devoted to address the highly variable and highly energetic emission from the 
blazars, while section [3731 presents results of a long-term study of gamma-ray burst afterglows. 

3.1. Parsec-scale jets 

Parsec-scale jets are radio features on a scale of several light years. They are seen emerging from 
AGNs over periods of months to years. It is widely accepted that the radio maps are not direct 
observations of the physical quantities (density, pressure, velocity, composition, etc.) in the jet, 
but are influenced by a number of relativistic effects (Doppler shift, beaming) as well as by the 
degradation of the image due to a finite resolution of the radio telescopes. By using numerical 
simulations of the jet dynamics and, subsequently, by computing the synthetic radio maps we 
can test theoretical jet models and directly compare them to the observations, as well as study 
some of the events which are expected to occur in thenjf]. In a previous work [9] we studied one 
such event, which is an injection of a perturbation into a steady jet shown on Figure [TJ This jet 
is over-pressured and under-dense with respect to the surrounding medium and has a number 
of recollimation shocks, which can be seen as knots in the density profile shown on Figure [TJ 

The perturbation has a form of a temporary increase in the fluid velocity at the jet injection 
point located at the left grid boundary. This results in the development of shock and rarefaction 
waves which propagate down the jet and interact with the recollimation shocks, resulting in a 
complex, non-linear evolution of the perturbation and several trailing components (see section 7 
of [9] or more details). Figure [2] shows the radio maps taken at six different epochs of the jet 
evolution. In the first three panels we can see the main component (perturbation) as a dark 
region propagating through the jet. In the last three panels we see a gradual reestablishment of 
the steady jet. 

The contours show the image degraded to the resolution available with today's radio 
telescopes. As can be appreciated, many of the small-scale details of the jet evolution and 
emission are lost in the relatively strong blurring of the image due to the low resolution of the 
actually observable radio maps. Therefore it is important to use the hydrodynamic simulations 
and possess the ability to compute time-dependent radio maps from those simulations because 
this enables us to compute both the full-resolution and the degraded images and thus evaluate 
which of the observed features of the jet dynamics and emission are intrinsic to the jet and which 
are the artefact of the finite observational resolution. 

3.2. Blazar jets 

Blazars are jets pointing almost directly towards us. They are characterised by high variability 
of their emission, often in the form of flares observed by X-ray satellite^. It is thought 
that the cause of the flares are the collisions of parts of the jet with different velocities 
("shells" in the following), whereby a fraction of the shell kinetic energy is dissipated by the 
shocks which propagate through the shells as a result of their collision, and a fraction of the 
dissipated energy is radiated and observed as a flare. This is called the internal shock scenario 
[EllISlEOlEIllSlEaEllSlElE]. Figures M and S illustrate the internal shocks model and the 
typical distance scales at which shells collide in blazar jets. 

We study the dynamics and emission from single shell collisions, using one-dimensional 
|23| and two-dimensional [21] hydrodynamic simulations, as well as one-dimensional 

4 Our work on blazar jets fsection l3.2[) as well as studies of the stability of initially very magnetized (a > 100) jets 
(e.g. |16] ) indicate that at parsec-scales distances AGN jets are expected to be at most very weakly magnetized 
(e.g. a < lO" 4 ) 

5 For an example of X-ray flares from a well known blazar Mrk 421 see e.g., [17] 
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Figure 1. Logarithm of the density of 
a steady-state jet propagating through a 
stratified atmosphere. The jet is initially 
under-dense and over-pressured with respect 
to the surrounding medium. The density is 
shown in arbitrary units, while the length 
scale is in the units of the jet radius at the 
left grid boundary. See section 2 of [9] for 
the details of the hydrodynamic models. Note 
that the vertical scale is enlarged for the 
better visualization of the jet. 
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Figure 2. Radio emission at 43 GHz re- 
sulting from a temporary velocity perturba- 
tion propagating through the steady-state jet 
shown on Figure [TJ The gray shades in the 
panels show the radiation intensity (in arbi- 
trary units) at (from top to bottom) 0.02, 
0.39, 0.75, 1.12, 1.94 and 4.58 years since the 
injection of the perturbation. The contours 
shows the image degraded to the resolution 
available with today's VLBI technique. The 
dashed blue line shows the approximate tra- 
jectory of the main component. See section 7 
of [9] for more details. 
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Figure 3. An illustration of the internal 
shocks model: an intermittently active central 
engine (supermassive black hole at the center 
of the galaxy) ejects inhomogeneous shells 
which collide and produce observed flares. 
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Figure 4. Distance scales in the internal 
shocks model in blazars: the shells, whose 
characteristic size is ~ 10 14 cm, collide at 
distances ~ 10 16 — 10 17 cm from the central 
engine, merge and continue flowing as a 
smooth jet at parsec-scale distances (« 10 18 
cm). 



magnetohydrodynamic simulations |SJ. We found that the lateral shell spreading is negligible 
on the length scales considered and that the one-dimensional approximation is accurate |21j . 



Furthermore, we find that the magnetization of the flow plays a crucial role in the efficiency 
of the conversion of kinetic into thermal energy (and subsequently into radiation). From the 
computational point of view we found the hydrodynamic simulations to be rather expensive and 
unsuitable for parameter space studies. 

In the follow-up works j3[ [13] we simplify the hydrodynamic approach in order to be able 
to devote more resources to the computation of the emission and to be able to cover larger 
parameter space. We study a large number of shell collisions where both shells have the same 
energy and size, but their velocity and the degree of magnetization can vary. We denote the shell 
Lorentz factor by T := (1 — v 2 /c 2 ) _1//2 , where v is the shell velocity. The magnetization for cold 
and relativistic shells (e.g. p/pc 2 <C 1 and r> 1, where p and p are the shell rest-mass density 
and the thermal pressure) is defined as a := B 2 / (AttT 2 pc 2 ) . We compute the exact solution 
of the Riemann probleiro for each shell collision and determine the strength and velocity of 
the shocks which are formed by the shell interaction. By covering a large parameter space of 
possible shell magnetization (we studied 10 6 shell collisions with a ranging from 10 -6 to 10 3 for 
each shell independently) we find that the "sweet spot" for the efficiency of the kinetic energy 
conversion is not for a <C 1, as might be expected (the smaller the a, the easier it is to shock the 
fluid), but rather for {gl = 1, 0\r = 0.1), where subscripts L and R denote the faster (left) and 
the slower (right) shell, respectively [3]. This result makes the magnetized internal shock model 
viable. In the follow-up work [13] we compute the multi-wavelength, time-dependent emission 
from ~ 10 3 shell collisions 
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Figure 5. Optical (5 x 10 14 Hz), X-ray 
(3 x 10 18 Hz) and 7-ray (10 21 Hz) flare (black, 
blue and red lines) produced by a collision of 
a faster shell with the initial Lorentz factor 
= 20 with a slower shell Tr = 10. 
The shell magnetization is &l = 10 -3 and 
an = 10 _1 . The shell geometry is cylindrical 
(radius 3x 10 14 m, length 6x 10 11 m) and their 
initial energy 10 52 erg. The flux is shown in 
units Jy Hz (1 Jy Hz = 10" 23 W / m 2 ). 



Figure 6. Average spectrum of the flare 
shown in figure [5j The total spectrum is 
shown in black, while the synchrotron, SSC 
and EIC contributions are shown in blue, red 
and green, respectively. The seed photons for 
the EIC are assumed to be monochromatic 
(frequency 2 x 10 14 Hz). The viewing angle 
(angle between the jet axis and the line of 
sight) is 5°, and the jet is assumed to be 
located at redshift z = 0.5. 



6 We use the exact Riemann solver [55] for the one-dimensional ideal relativistic magnetohydrodynamic Riemann 
problem with magnetic fields perpendicular to the flow velocity. 

7 We could not compute the emission from all of the 10 6 shell collisions studied in [3J because it was not feasible 
with present supercomputers (the computation time for the emission from 10 shell collisions exceeds 200 thousand 
hours) . 



Figures [5] and [6] show an example of the light curve and spectra produced by a collision 
of two moderately magnetized shells in a typical blazar jet (Rueda, Mimica & Aloy, 2012, in 
preparation). In figure [5] can see that the optical and 7-ray light curves peak somewhat earlier 
than the X-ray light curve. The reason is that the SSC emission (section I2.2.2P is produced by 
the scattering of the synchrotron emission which needs a finite time to travel from the place 
where it is produced to the site of its scattering. The sharp drop in the emission is related 
to the moment when shocks cross the shells and cease dissipating energy. Average spectra are 
shown in figure EJ The low- frequency emission is dominated by the synchrotron emission, while 
the high-energy spectrum is dominated by the SSC component. As can be seen, the correct 
computation of the SSC component is essential, and since this component is computationally 
most expensive, it puts an upper limit on a number of models which can be computed on today's 
machines (see footnote [7|) . 

3.3. Gamma-ray burst afterglows 

Gamma-ray bursts (GRBs) are produced by an internal dissipation of energy in the relativistic 
outflow. The GRB afterglow emission is produced by the interaction of the relativistic jet with 
the external (circumburst) medium. It is not known to what extent the jet is magnetized. If the 
outflow is initially dominated by the thermal energy (fireball model |26t 127] ) then it is expected 
that its a <C 1. If, on the other hand, the outflow is initially dominated by the Poynting-flux 
[28} 129] . then it is possible to find a significantly magnetized jet at the onset of the afterglow 
phase. 

We study the interaction of an arbitrarily magnetized GRB jet with the external medium. 
We model the interaction by considering the deceleration of a dense, magnetized shell as it moves 
through and compresses a constant density medium in front of it. One of the crucial differences 
between the evolution of a non-magnetized and a magnetized shell is in the existence of a reverse 
shock (a shock propagating through the shell): if a is high enough, then the interaction between 
the shell and the external medium is not strong enough to shock the shell material itself. By 
analyzing timescales of the shell deceleration and the propagation of magnetohydrodynamic 
waves through the shell we find the conditions for the existence of a reverse shock depending 
on the GRB magnetization, energy, Lorentz factor, duration and the density of the external 
medium [30]. Figure [7] shows the parameter space spanned by the ejecta magnetization and a 
dimensionless parameter £ (see figure caption for definition). As we can see, there is a substantial 
portion of the parameter space (shaded region), where the reverse shock does not form. Since 
it is absent, the particles are not accelerated and we expect to fail to observe an emission 
from this shock during the early afterglow, which would be consistent with the absence of such 
observations in most GRB optical afterglows |31j . 

In the follow-up works [32] we confirmed these results by means of numerical simulations. 
Using MRGENESIS we performed high-resolution, long-term one-dimensional relativistic 
magnetohydrodynamic simulation of a shell-medium interaction, and then use SPEV to compute 
the resulting light curves. Figure [8] shows an example of the optical light curve produced by non- 
magnetized, weakly-magnetized and strongly-magnetizeed shells with the same initial energy, 
Lorentz factor and width. As can be seen, the strongly magnetized shell light curve does not 
possess the optical flash feature that the other two possess. Therefore, strong magnetization 
might be required for the paucity of the observed optical flashes to be explained. Since most of 
the observed GRBs have £ in the range [0.3, 3], this would require a > 0.1 for almost all GRBs 
[301 [32]. This would imply that GRB jets could be strongly magnetized even at such late stages 
of their evolution as the afterglow phase, which puts constraints on the amount of magnetic 
dissipation that can happen at previous stage (prompt emission). 

The simulations needed to correctly compute afterglow light curves are very computationally 
demanding due to the high resolution required in the direction of the jet motion (e.g. more 
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Figure 7. Parameter space of 

the ejecta magnetization <to and 
t := (3^/4vrp ext c 2 ) 1/6 Ao 1/2 r o 4/3 , where E, 
Ao and Tq are the ejecta initial energy, width 
and the Lorentz factor, and p ext is the ext. 
medium density. The shaded region shows 
the ejecta which do not form a reverse shock, 
while the vertical boundary at ao = 0.3 
indicates those ejecta whose radial spreading 
becomes significant and makes the formation 
of a reverse shock inevitable 1 30 1 . 



Figure 8. Optical (5 x 10 14 Hz) light curve 
produced by the non-magnetized (<7q = 0), 
weakly magnetized (do = 0.01 and strongly- 
magnetized ((To = 1) shells (black, blue 
and green lines, respectively) whose initial 
conditions are E = 3 x 10 47 J, T = 300 
and Ao = 3 x 10 9 m. The non/weakly 
magnetized ejecta have a strong reverse shock 
which produces an optical flash, while the 
strongly magnetized ejecta light curve does 
not have this feature |32j. 



than 10 4 zones in the initial shell) and the long timescales needed to compute the light curve 
of sufficient duration (approximately 10 7 numerical iterations). Therefore, performing this type 
of calculations to cover a full parameter space of feasible GRB afterglows is not possible at 
the moment. However, the existence of re-scaling laws (see section 4.4 of [7]) can help us to 
perform a small number of simulations and then extend their results to cover a somewhat larger 
parameter space. 

4. Conclusions 

We have presented a framework consisting of a numerical relativistic magnetohydrodynamic code 
MRGENESIS and the non-thermal radiative transfer code SPEV. We show how its modular 
nature has enabled it to be successfully been applied to a number of current problems and 
issues in the relativistic jets physics, and how it enables us, via the computation of synthetic 
observations, to directly compare the results of our simulations with the observations from 
ground- and space-based telescopes. In the future we plan to improve the scaling of MRGENESIS 
to more than 10 4 cores, and also to add new radiative processes to SPEV (e.g., polarization, 
improved inverse-Compton, radiative transfer in curved space-times, etc.). The code scaling 
improvements will allow us to perform full 3D magnetohydrodynamic simulations and radiative 
transfer calculations, relevant for astrophysical scenarios where axial symmetry cannot be 
assumed. 
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